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Abstract 

The parameters of many-body potentials for Co, Nb and Zr metals, based on the embedded- 
atom method, have been systematically derived. The analytical potential scheme allows us to 
reproduce correctly the cohesive energies and structural properties of the pure metals and selected 
alloys making use of a small set of parameters. With a pair potential going smoothly to zero 
for a sufficient cutoff radius, radial partial and bond angular distribution functions for Co, Nb, 
Zr and alloys are computed using molecular dynamics simulations that ensure good quantitative 
agreement with the available experimental data up to the melting point. Atomic short range order 
is analysed in the light of consecutive Gaussian function decomposition and Honeycutt-Andersen 
indices. 

PACS numbers: 61.25.Mv, 71.15.Mb 
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I. INTRODUCTION 



Today we are not surprised that a non- crystalline solid orders magnetically. It is known 
that, with few important exceptions, the amorphous and crystalline phases of the same 
material do not differ very much magnetically. In the last three decades, the discovery 
and large scale investigation of rapidly solidified alloys have made possible a new scenario 
in basic and applied magnetism. The unusual behaviour of the bulk metallic magnetic 
glasses (BMMGs) - formed by supercooling the liquid state of certain metallic magnetic 
alloys - usually occurs for systems containing atoms that exhibit a well-known sensitivity 
to the immediate neighbourhood . However, despite the large number of technological 
applications of BMMGs, the detailed origin of the links between structural and magnetic 
properties has yet to be established, but is strongly dependent on the atomic short-range 
order j^. 

Nowadays, atomic-scale simulations of solids based on interatomic potentials are routinely 
performed to explore the short-range order of BMMGs 0, [sl as a precursor to the study of 
the magnetic ordering, which is evidently beyond such methods. However, the construction 
of realistic n-body potentials is mandatory to any simulations. When applied to metals and 
alloys, several major methods and extensions of construction of such potentials have been 
established from density functional theory, i.e., the embedded-atom method (EAM) 6| or 
tight-binding second moment approximation, i.e., the Finnis- Sinclair (FS) model 7| and 
related models j^, . The two models have very similar computational requirements, and 
the names are often used interchangably, however there are some distinctions which come to 
the fore when considering multicomponent alloys [10]. Because of the parameters involved 
in these model potentials, the EAM method was first used to study simple metals in their 
relevant crystalline structur es [ill . Il2l|. Extensions to binary alloys were formulated with 
special attention to hcp-fcc ij] and hcp-bcc 15| systems. Because they are often based on 
spin-less approximations of electronic density, such effective model potentials traditionally 



n 



neglect the magnetic ordering Il6l . Il7l |. Notable exceptions have been published 



recently by Dudarev et al. [l8| and Ackland et al. 19|, |20], but such methods are still in their 
infancy. Even in spin-less schemes, these effective potentials may be able to represent the 
compositional ordering which is, at least in the localised magnetism picture, a prerequisite 
for understanding the magnetic behaviour. Existing magnetic potentials concentrate on the 
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one-site magnetism, for which the energetics reduces in form to a simple embedded-atom- 
type potential, explaining the success of standard EAM schemes on magnetic materials. 

Supercooled Coi-^^Nhj; and Coi-ajZrj, magnetic metallic glasses in the cobalt-rich region 
represent test cases for both hcp-bcc and hcp-hcp metal sub-systems. These compounds 
have been experimentally studied because they are strong ferromagnets, as revealed by a 



high value of the exchange constant [2l|. Consequently, large short-range compositional 
inhomogeneities should induce significant variations of the long-range magnetic order. The 
present work considers both the construction of embedded-atom potentials for such materi- 
als, as well as the application to molecular dynamics (MD) simulations in order to assess one 
of the issues in the large field of these BMMGs, namely short-range order quasicrystallinity. 

II. METHODOLOGY FOR POTENTIALS 

In the FS method, the total internal energy E of a A^-atom system and the electron 
density, p(Ri), for an atom located at Rj due to all other atoms are given as; 

^ N N 

N 

p(R.) = X;/M, (2) 

where f{rij) is the electron density at atom i due to atom j as a function of the distance 
between them, = ||Rj — Rj|| is the separation distance between atoms i and j, F(p(Rj)) is 
the energy to embed atom i in an electron density p(Rj), and 4>{rij) is a two-body potential 
between atoms i and j. As long as an angular independent formulation is considered, the 
electron density is a radial function only. For an alloy model, an embedding function, F, has 
to be specified for each atomic species supplemented by an atomic electron-density function, 
/, and a two-body potential, 0, specified for each possible combination of atomic species. 
For uncompressed metals, Gupta [3] and Tomanek et al. [23^ have shown that the host 
electron density can be represented as an exponentially decreasing function of the distance 
to better account for atomic relaxation near impurities and surfaces. In this approximation, 
/ is given as; 

f{r) = feexp{-x{r/re-l)), (3) 
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where /e is a scaling factor determined by the cohesive energy, Ec, and the atomic volume, 
Te is the nearest-neighbour distance in the relevant pair of atoms and x is an adjustable 
parameter. Analysing the interatomic interactions in effective-medium theory, Jacobsen et 



al. 2J] have shown that if an exponential form is chosen for the density function, then the 
interatomic potential, 0, should also be an exponential function of the distance. In this 
study, the interatomic potentials of all the pairs considered are defined by a potential of the 
form; 

74exp(— r/ro) < r < ri, 

5 

^a^r* n < r < r^, (4) 

< r, 

where the interaction is designed so as to go smoothly to zero at the distance according 
to a polynomial spline function. The potentials are constructed subject to the constraints 
that the radial functions and their first and second derivatives must be continuous at the 
boundary points, and also that the function must have a stationary point at r^- Once 
A, ro, ri and are fixed, this procedure ensures that the coefficients, {ctj}, are uniquely 
determined by solving a simple 6x6 linear system of equations. These coefficients are 
reported in Table [T] for completeness. For all the pairs of atoms, = 2.5A is kept fixed 
and cor resp onds to a typical radius where the stiff repulsive part of the potential ends in 
metals 



25 



261, 



271]. Because of the screening in metals, the stationary point is located at 
least between the second and third nearest neighbours for the lowest energy crystal phases 
as previously noted [28]. Since the embedding energy is assumed to be independent of the 
source of the electron density and the hopping integrals are a function only of a radial 
distance between atoms, the embedding functional F is taken as 

F[pir)] = vW)- (5) 

This functional form gives a band energy proportional to the square root of the second 
moment of the electron density of states [22]. However, moments of higher order cannot 
3e expressed in such a simple analytical form and a more complex method must be applied 



291 ]. Johnson 13| has considered that since the electron density at any location is taken 
as a linear superposition of atomic electron densities, this function should be taken directly 
from monoatomic models with a relative scaling factor between elements for an alloy model. 
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On the other hand, Finnis and Sinclair [3] and Cleri and Rosato [s^ have considered mixed 
pair electron-density functions not necessarily connected to the atomic ones, removing the 
alloy scaling factors. 

For hcp-Co and hcp-Zr, the parameters for the atomic electron-density and the inter- 
atomic potential are fitted in order to reproduce the experimental cohesive energies, the 
unit cell parameters and the five independent elastic constants of these systems as given by 
Cleri and Rosato [s^. Moreover, cell parameters and elastic constants of fcc-Co and bcc- 
Zr are also included during the fitting procedure as taken from references 



32 



and 



references therein. Generally, for a small cutoff distance the largest number of interacting 
neig hbours per atom in the crystalline structure leads to the more stable phase. Ducastelle 



3J] has shown that in the second-moment approximation, with interactions restricted to the 



nearest neighbours, the cohesive energy for the hep and fee phases is the same and the c/a 
ratio is equal to the ideal value. Hence, it is necessary to go up to at least the fourth-moment 
approximation to discriminate between the fee and hep phases and to give a value of c/a 
different from 2 a/2/3. In our case, since the c/a ratios are not taken to be the ideal one and 
because the potentials have a very short range, the cutoff distance of atomic electron density 
should be larger than that of the potential. So for the electron density, the cutoff distance 
is taken to be 4.87A, thus including up to seven shells of neighbours within hcp-Co, and 
three for both hcp-Zr and bcc-Nb, which are all the stable phases for each pure element. A 
check is also performed to ensure that at least three shells of neighbours have been included 
for the high pressure/high temperature structures since this is necessary to keep the relative 
energies of each phases in the correct order 33|. Moreover, it has been observed that in 
incorporating the elastic constants of the bcc-Zr phase in the database, the fitting of both 
the cell parameters and cohesive energy of the hcp-Zr is rather poor with this cutoff radius, 
so only the hcp-Co, fcc-Co, hcp-Zr and bcc-Nb elastic constants are included during the 
fit. For Nb, the cohesive energy, the lattice parameter and the three independent elastic 
constants of the bcc-Nb phase are taken from reference 35|] and the theoretical fcc-Nb cell 
parameter and cohesive energy are also included 36|]. No elastic constants of the fcc-Nb 
were found to incorporate into the training set of observables. 

The selection of the functional form taken in Eg. ([!])- ([2D i s extended to AB alloys based 



on the second-moment form that has been applied to Zr 



37|. The embedding function. 



atomic electron-density function and two-body potential are assumed to be of the same 
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form as in Eq.([3])-(jl]), with i^ab and (pBA assumed to be equal. The alloy potentials and 
atomic electron-density functions are determined independently of the monoatomic counter- 
parts if sufficient data are available. However, it is known that for equilibrium immiscible 
systems it is a challenging task to fit cross potentials, since there is often insufficient ex- 
perimental data related to the respective alloy compounds. In order to circumnavigate 
this problem, density functional calculations have been per formed on selected intermetal- 
lic structures using the Quantum-ESPRESSO package |38|]. For these calculations, non- 
local ultrasoft pseudo-potentials are employed in combination with a plane-wave basis set. 
The generalized-gradient approximation, as parametrized by Perdew, Burke and Ernzerhof 



39( 1 . is selected for the exchange and correlation term. For the Brillouin zone sampling, a 
12x12x12 Monkhorst-Pack mesh is used for the fc-point summation in the self-consistent 
calculations [40] for all the primitive cells, which leads to converged structural parameters 
to within 0.1% of the cell parameter. Thus the lattice constants and cohesive energies of 
several Co-(Zr,Nb) and Zr-Nb crystalline structures reported in Table ( IIIII) are obtained and 
then included in the fitting procedure for the Co-(Zr,Nb) and Zr-Nb cross potentials and 
atomic electron-density functions. The parameters of the fitted terms are listed in Table [B 
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') 


-0.5668 


-1.6874 
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-6.3808 



TABLE I: Potential and atomic electron-density parameters for (Co,Nb,Zr) systems. 
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The overall fitting procedure is performed in two separate steps. First, the densities 
and potentials are derived for the simple metals. Once the corresponding parameters are 
obtained, the fit is applied to selected binary metals for cross-densities and potentials without 
altering the terms for the simple metals. The same cutoff radius is kept constant during all 



the steps and the calculations are performed within the GULP computer code 41 1- 

In Table [Tll a list of some basic physical properties as computed by the present set 
of potentials and the corresponding experimental values are shown for Co, Nb and Zr. 
The fit correctly reproduces the structures and properties of hcp-Co, fcc-Co, hcp-Zr and 
bcc-Nb. The absolute average percentage difference between calculation and experiment is 
found to be 0.6%, 0.7% and 8%, for the cell parameters, cohesive energies elastic constants, 
respectively. However, the elastic constants of fcc-Nb produced a negative Young's modulus, 
as expected, in this excited locally unstable structural phase 4J]. The properties of bcc- 



Zr are reproduced with a sufficient accuracy, except for the C44 elastic constant where the 
largest percentage error of 45% occurs. The potentials fitted by Willaime et al. [33| also 
exhibit such a discrepancy though with a much larger error. This may be corrected by 
relaxing the constraint of the square-root form of the embedding functional and considering 



more neighbours 12| or by including explicit angularly dependent terms in the potentials 
at a cost of additional parameters [45| . 

In Table llllt a list of some basic physical properties fitted from these potentials and 
the corresponding ab initio calculated values for selected binary alloys is shown. The B2, 
C15 and LI2 denominations are relative to the strukturbericht structural types classification 



481], such as B2 is the CsCl structure type, C15 is the MgCu2 structure type and LI2 is 
the CusAu type. The fitting procedure appears to correctly reproduce the ah initio derived 
phase stability order and lattice parameters. However, the absolute cohesive energies are 
reproduced to a lesser extent. This may be also improved consistently by increasing the 
cutoff radius on Zr and Nb electron density terms, which includes more neighbours in the 
total energy sums of Eq.([l]). 

III. APPLICATION TO SIMPLE METALS 

The validity of this potential in describing the atomic interactions can be illustrated 
outside the original systems used for parametrisation by considering liquid cobalt. An MD 
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TABLE II: Physical properties for Co, Zr and Nb simple metals as fitted with a cutoff radius of 
4.87A. SI) Cohesive energies and lattice parameters are taken from Kittel S-^) elastic constants 



are taken from Simmons and Wang 



43|. 



simulation with 300 atoms, which allows an individual description up to the 13* nearest- 
neighbors in hcp-Co, and periodic boundary conditions was first performed within the iso- 
baric, isothermal ensemble (NPT) at a temperature of 1670 K, which is slightly below the 
experimental melting point. The simulation was run for 200 ps with a time step of 0.1 fs to 
simulate the radial pair distribution function (RPDF). Once thermal equilibrium is reached, 
the RPDFs are sampled every 0.2 ps to produce an average. As these RPDFs are subject 
to statistical noise, a smoothing formula is used to replace each RPDF by a least-squares 
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TABLE III: Lattice parameters and cohesive energies for selected binary metals in their corre- 
sponding symmetry as computed with a cutoff radius of 4.87A. (1) experimental value of 3. 181 A 



46 1, (■^) experimental value of 6.774A 



47[. T 



le B2, C15 and LI2 classification is relative to the 



strukturbericht structural types classification 



48] 



polynomial that fits a sub-range of several points. For all the calculated radial functions, 
a third-degree, five-point smoothing procedure is applied several times on the data until 
convergence 49]. MD simulations are repeated on three different atomic configurations at 
the same temperature in order to sample more accurately the configurational space and 
an overall averaged RPDF is computed. This averaged RPDF compares very well to the 



experimental data 



50| as shown in Fig. [H The agreement is similar to that presented by 



Bhuiyan et al. 5l| and more recently by Han and coworkers 52] using much more elaborate 
EAM potentials, but in these works no quantitative analysis of the RPDFs was performed. 
However, the first peak of the RPDF g{r) appears to be asymmetric and might be com- 
posed of more than one atomic shell. Following Kita et al. 53(], a decomposition of rg{r) 
in consecutive Gaussian functions is applied. The insert in Fig. [1] shows the decomposition 
of rg{r) at 1670 K in six Gaussian functions where all 18 parameters are allowed to vary 
freely during the fit. The first coordination shell is defined by a cutoff distance r^, which is 
taken to be the first minimum of g{r). For this temperature, = 3.46A. The first peak is 
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composed of three Gaussian subpeaks located at ri = 2.407A, r2 = 2.645A and = 3.095A, 
respectively. The coordination number A^^ is calculated by integrating the Gaussian function 
according to A^^^ = ATr^^'^UQAari, where uq is the atomic density, A is the amplitude of the 
Gaussian function, a is the square root of the variance, and is the maximum radius. For 
each subpeak, Nc is equal to 3.46, 4.71 and 3.81 with a sum of 11.98. For a temperature of 
1800K, which is slightly greater than the experimental melting point, the sum decreases to 
11.57. These coordination numbers are close to the experimental ones found in liquid Co [s^ 
(12.5±0.5 at 1670 K and 12.1 ±0.5 at 1800 K) and consistent with those calculated for other 
metallic systems 5J|. Such a high value of the coordination number and the possibility of 
decomposing the first peak of g{r) is an indication that the short-range order of the liquid Co 
is more complex than the one given by a simple icosahedral ordering as suggested by Holland 
et al. js^. Moreover, performing MD simulations for these two temperatures allows us to 
predict a variation of the density with the temperature of dp/dT = —9.68 10~^gcm~^K~^ in 
good agreement with the experimentally reported value of dp/dT = —9.88 10~^gcm~^K~^ 
value [48 1 . 

MD simulations were repeated under the same conditions for pure Zr to a higher tem- 
perature of 2290 K, above the experimental melting point. The experimental RPDF is 
compared against the simulation of this liquid state, shown in Fig. [2l For this metal at that 
temperature, the first peak is composed of two Gaussian functions located at ri = 3.09lA 
and r2 = 3.574A, respectively. The ratio r2/ri = 1.156 is close to that of the two first 
nearest-neighbour distances for a bcc lattice 2/^/3 = 1.1547. The coordination number 



for each subpeak is 6.10 and 5.59 wit 
one of 11.9 ± 0.5 found in liquid Zr 



1 a sum of 11.69. This value is close to the experimental 



551]. Even if the ratio of the radii tends to favour a bcc 



lattice, the corresponding coordination is very different. 

MD simulations are repeated under the same conditions for pure Nb to a higher tem- 
perature of 2750 K, which is the experimental melting point. The simulated RPDF of this 
liquid state is shown in Fig. [3l For this metal, the first peak of rg{r) is strongly asymmet- 
ric and is composed of three Gaussian functions up to Vc = 4.0A, located at ri = 2.69lA, 
r2 = 3.14lA and = 3.840A. The ratio r2/ri = 1.167 is greater than that of the two first 
nearest-neighbour distances for a bcc lattice. The coordination number Nc for each subpeak 
is 2.95, 7.06 and 3.10 with a sum of 13.11. To our knowledge, no experimental data on the 
radial distribution function of pure Nb is available in order to compare with. 
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The short-range order can also be examined by calculating the bond angle distribution 
functions g{6) that represent the angle between the bonds connecting a central atom to 
two neighbouring atoms, as illustrated in Fig. IHfor liquid cobalt, zirconium and niobium. 
The angle is calculated for pairs of interatomic distances given by a cutoff corresponding to 
the first minimum of the RPDF (i.e. 3.5A for Co, 4.3A for Zr and 4.lA for Nb). In the 
case of Zr, the calculated distribution exhibits a prominent peak near 6 = 57° (close to an 
equilateral triangle), a broader maximum near 6 = 109° and a rather flat maximum near 
6 = 150°. In the case of Co, the first peak is broader and close to 55° whereas the second 
peak enlarges but remains at 109°. For Nb, the situation is different with two broader peaks 
near 50° and 99°. Viewing the structure in terms of dominant clusters, the bond angle with 
the highest density at the nearest neighbour distance are for a regular icosahedron 63.4° and 
116.4°, while for fee the prominent angles are 60°, 90° and 120°. For hep, angles of 109.471° 
and 146.443° are added when the ratio c/a = 2 a/2/3 but they are less frequent. For a bcc 
lattice, the prominent angles are 70.53° and 109.471°. In the Zr structure, the first peak 
tends to favour the fee and hep structures while the second peak tends to favour the hep 
and bcc. This means that the dominant structure should be hep. However, the angles of 
90° and 120° are not so strong whereas some defective icosahedron angles should be there 
too. This suggests a predominantly distorted icosahedral character. For the same reasons, 
the case of Co also favours defective icosahedron as well. For Nb as the first peak is much 
more located near 50°, this suggests much more intriguing short-range structures with less 
neighbours. 

To assess more quantitatively the local structures in amorphous alloys, Honeycutt and 
Andersen (HA) analysis has been proven to successfully differentiate face-centered cubic, 
hexagonal close packed, icosahedron and binary bcc structures [s?]. To perform such 
analysis, a set of four indices is constructed for each pair: (i) the first index denotes to 
what peak of the RPDF, g{r), the pair under consideration belongs; (ii) the second index 
represents the number of near neighbours shared by this pair; (iii) the third index counts the 
number of nearest-neigbour bonds among the shared neighbours; (iv) and a fourth index is 
used to differentiate configurations with the same three indices, but with a different topology. 
For instance, fee crystals are fully described by four pairs, such as 1421, 2101, 2211 and 2441, 
whereas hep crystals also contains the 1422 and 2331 pairs in addition. Moreover, the 1441, 
1661, 2101, 2211, and 2441 are the only pairs in a perfect bcc crystal. Icosahedral order is 
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described by Mackay icosahedra, composed of twinned, distorted fee tetrahedra with an index 
of 1551, whereas the 1541 and 1431 indices are more characteristic of a distorted icosahedral 
local order. Up to the distance cutoff corresponding to first minimum of each g{r), Table HVl 
reports the HA analysis in the liquid state of each simple elements. The microscopic analysis 

index Co Zr Nb 

1311 0.07±0.01 0.06±0.01 0.03±0.01 

1321 0.07±0.01 0.06±0.01 0.05±0.01 

1421 0.04±0.01 0.03±0.01 0.02±0.01 

1422 0.07±0.01 0.06±0.01 0.04±0.01 
1431 0.19±0.01 O.lTibO.Ol 0.13±0.01 
1441 0.04±0.01 O.OSibO.Ol 0.07±0.01 
1541 0.15±0.01 O.lSibO.Ol O.llibO.Ol 
1551 0.12±0.01 O.llibO.Ol 0.14±0.01 
1661 O.OSibO.Ol 0.06ib0.01 0.07ib0.01 
2101 1.58ib0.02 1.55ib0.01 1.52ib0.02 
2211 0.95ib0.02 0.94ib0.02 0.96ib0.02 
2321 0.23ib0.01 0.22ib0.01 0.21ib0.01 
2331 0.50ib0.02 0.53ib0.02 0.58ib0.02 
2441 O.OSibO.Ol O.OSibO.Ol 0.07ib0.01 

TABLE IV: Honeycutt and Andersen analysis of the simulations in the liquid state for Co 
(T=1670K), Zr (T=2290K) and Nb (T=2750K). 

emerging from the data of Table [IV] indicates that the short-range order of the liquid state 
is dominated by distorted icosaheral and icosahedral structures since the 1541, 1431 and 
1551 indices respectively are large as anticipated. The high value of the 2331 pairs is also 
an indication of the icosahedral order. The small distortion from perfect icosahedral order 
observed in the angular distributions of these liquids suggests that the local icosahedral 
order should dominate. However, small distortion form a perfect tetrahedron does not form 
different HA indices from those for a perfect icosahedron. Our HA analysis shows that the 
icosahedral distortion is larger than reported by the bond angle distribution curves. This 
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result is in agreement with the experimental investigation on liquid Ti, Zr and Ni conducted 



by Kim and Kelton 
Pasturel 



581]. Using first-principles molecular dynamics simulations, Jakse and 
59( 1 have concluded there exists competition between a polyhedral and bcc-type 
short-range order in liquid and supercooled Zr, whereas Kim and Kelton 3] have reported 
no regular dominant cluster type that can describe the experimental liquid structure of 
transition metals, including Zr. This is supported by the values of the HA indices reported in 
Table HVl which are not very different from each element in the liquid state. The abundance 
of the 1661 pairs indicates that bcc order is very low, but slightly increases when going 
from Co to Nb. Interestingly, the lowest energy geometrical structures of magnetic cobalt 
clusters mainly follow an icosahedral growth pattern with some cubic-type structures at 



some particular sizes 



60|. 



For liquid Zr, these low values have been re por ted both experimentally 
first-principles molecular dynamics simulations 



and using 



59|. Furthermore for liquid Zr, Jakse and 



Pasturel have performed HA analysis on the inherent structures and found an abundance of 
the 1551 pairs in the liquid state that is twice the value found here. Such structures indicate 
the presence of perfect icosahedra as a local minima of the potential energy surface. For 
instance on liquid Nb, the 1551 index goes to 0.30±0.01 on inherent structure. 

IV. APPLICATION TO BINARY ALLOYS 

To validate the quality and transferability of our potentials, lattice constants and atomic 



cu. 



61 



ated 



62|. 



internal positions of several crystal structures not entering in the fit have been ca 
minimising the free energy at T=300K, and compared with the experimental values 
The results are summarised in Table |V] for the varying degrees of freedom according to the 
corresponding space group and the agreement is generally good with an average absolute 
error of 1.76%. To complete the validation, the elastic constants of CoNb and CosNb 
are calculated ab initio applying finite differences to the stress tensor and compared with 
those calculated analytically by our potentials. The results are shown in Table |Vl] and the 
agreement is satisfying. As for liquids, MD simulations are performed on Coo.gZro.i for a 
300 atom system with periodic boundary conditions. First, the atoms are placed randomly 
into the simulation cell using a hard sphere criterion based on their atomic radii and the 
cell volume is then adjusted according to the phenomenological Miedema theory [63], which 
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CoyNbg Co23Zre 



exp.[61j EAM(300K) exp.[62] EAM(300K) 

a(A) 5.01 5.019(+0.18%) a(A) 11.516 11.484(-0.28%) 
24.66(-6.95%) xcog 0.378 0.3791(+0.30%) 
0.5005(+0.10%) XC704 0.178 0.1732(-2.71%) 
0.4995(-0.10%) yizr^ 0.208 0.2075(-0.26%) 
0.5801(-1.67%) 



c(A) 26.5 
xcoi 0.5 
ycoi 0.5 
zcoi 0.59 



ZNhi 0.167 0.1664(-0.34%) 
7.Nb2 0.346 0.3211(-7.19%) 
Z7Vb3 0.448 0.4359(-2.71%) 



TABLE V: Lattice constants and atomic internal positions of CoyNb^ an d Co23Zr6 calculated with 
our potentials at 300K and compared to the experimental values 



61 



63]. 



CoNb {D2) CosNb (L12) 
ah initio EAM ah initio EAM 
Cii (GPa) 251 242 368 357 
C12 (GPa) 173 143 164 194 
C44 (GPa) 60 71 160 131 

TABLE VL ab-initio elastic constants of selected binaries as compared to the present set of poten- 
tials. 



gives good estimates for the experimental volume of glasses of these alloys. Then at constant 
pressure and a temperature of 1800 K (higher than the liquidus phase boundary for this alloy 
composition), MD has been run for 200 ps to simulate the liquid phase. The sample is then 
quenched at a rate of 7.5 x lO^^Ks"^ and maintained at 300 K for at least a further 200 ps. 
In Fig. [5l the simulated RPDFs are computed and compared against experiment 6J]. The 
simulated partial distribution functions for cobalt compare well to experiment including the 
position of the first and second peaks. However, the experimental Zr-Zr RPDF does not 
exhibit a structural trend, whereas our simulation does. RoBler and Teichler have reported 
similar results in their study of atomic mobilities and structural properties of supercooled 
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amorphous Coi-^Zr^. using very different interatomic potentials {4]. This suggests a possible 
lack of resolution both in X-ray diffraction and EXAFS because of a low contribution in the 
spectra of these minority atoms, as previously anticipated 4|]. However, our potentials seem 
to reproduced correctly the position of the first peak for Zr-Zr and the double peak character 
of the second peak in the Co-Co and Co-Zr distributions in comparison with ref. J]. The 
maximum of the Co-Co peak is simulated to occur at 2.44A, compared with an experimental 
result of 2.42A [64]. Using the Gaussian subpeak analysis, this first peak is found to be 
composed of 3 shells at 2.43A, 2.55A and 2.88A with coordination numbers of 4.42, 4.83 
and 2.30. The total coordination number up to Vc = 3.12A is equal to 11.55 in comparison 
with 10.90 found experimentally. In our simulations a narrower scattering of the positions of 
the first 3 subpeaks is observed, in contrast to the liquid state, which is anticipated during 
the cooling. In the closest crystalline form, Co23Zr6, the highest coordination number is 
obtained for a pair of cobalt atoms located at 2.4334A in a cubic cluster and other local 
structures with 4 and 3 neighbours are also found at 2.36A, 2.5lA, 2.5lA and 2.8lA 62|. This 
suggests tetrahedral clustering below Vc or defective icosahedra up to Tc in this supercooled 
alloy. On the other hand, the maximum of the first peak in Co-Zr is simulated to be at 
2.79A, in excellent agreement with the reported experimental distance of 2.79A 64|. The 
analysis through Gaussian functions reveals 2 subpeaks located at 2.7lA and 3.00A with a 
coordination number of 8.68 and 8.40, respectively. This also suggests a clustering of the 
bcc-type at very short range. 

MD simulations have also been performed for Coo.gNbo.i under the same conditions. The 
reported RPDFs are shown in Fig. [6l In Coo.gNbo.i and Coo.gZro.i, the Co-Co first peak 
distance is calculated to be at 2.437A and is not affected by the non- magnetic added atoms 
at such a low concentration. The situation changes for the second and third peaks with 
a much more distinct third shell in Coo.gNbo.i than in Coo.gZro.i. Interestingly, the Co-Nb 
(resp. Co-Zr) first neighbour equilibrium distance is 2.49A (resp.2.79A). This is lower than 
the simple prediction related to their corresponding atomic radii (2.70A (resp.2.85A)) 65|. 



However, Jamet et al. 



661] have reported a Co-Nb distance of 2.58A in studying cobalt 



nanoparticles embedded in a niobium matrix. The simulated correlation in the minority 
pairs is structured in the Co-Nb system, as for Co-Zr. However, in both cases the first shell 
of neighbours seems to be depleted of their atoms to fill the second or third shells. It is 
doubtful that these structures are an artefact of the low number of atoms considered in our 
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simulations because RoBler and Teichler have simulated systems more than twice the size 
of ours and found the same behaviour. Up to Vc = 3.1A, the bond angular distributions 
of these two systems are calculated and shown in Fig. [71 These distributions exhibit well- 
structured peaks suggesting more crystalline environments including a 150° distinct angle 
in the Coo.gNbo.i- As the concentration of minority atoms is low, these distributions are 
dominated by hcp-like Co clusters upon cooling. This trend is more pronounced for added 
Nb atoms than Zr atoms probably because of their corresponding atomic radii. To assess 
such hypothesis, the HA indices are calculated and reported in Table IVlII In Coo.gZro.i the 



index C 


oo.gZro.i 


Coo.gNbo.i 


1311 


0.02 


0.06 


1321 


0.03 


0.03 


1421 


0.02 


0.06 


1422 


0.04 


0.11 


1431 


0.14 


0.21 


1441 


0.05 


0.01 


1541 


0.15 


0.17 


1551 


0.24 


0.16 


1661 


0.08 


0.03 


2101 


1.49 


1.57 


2211 


0.81 


0.87 


2321 


0.20 


0.20 


2331 


0.70 


0.62 


2441 


0.10 


0.13 



TABLE VII: Honeycutt and Andersen analysis of the simulations for supercooled CoQ.gZro.i and 
Coo.gNbo.i at 300K. The absolute error bars of the abundances are 0.01. 

dominant pairs are 1551 and 1541 exhibiting more icosahedral than distorted icosahedral 
order. This tendency is inverted in Coo.gNbo.i system with a larger amount of 1431 pairs 
favouring distorted icosahedral order. This is explained by a smaller difference in atomic 
radii between Co and Nb than Co and Zr because an atomic size difference of approximately 
10% can relieve spatial frustration and stabilise the icosahedral structure [67,] • In Coo.gNbg.i 
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the 1422 pairs are abundant indicating an hep order which is less present in Coo.gZro.i system. 
Even if the value is rather low, the 1661 pairs are of importance indicating a slight bcc order, 
as anticipated. It is observed that while useful for gaining understanding of the evolution of 

the dominant short-range order, the single cluster model cannot capture the richness of the 
supercooled binary alloys structures. 

V. CONCLUSION 

A fitting procedure has been performed to consistently derive a self-consistent set of many- 
body parameters for Co, Zr and Nb simple metals and selected alloys, including validation 
against first principles results where there are gaps in the experimental data. Combined 
with MD simulations, these parameters allow us to calculate RPDFs and bond angular 
distributions in the liquid phase for Co, Zr and Nb. Applied to supercooled binary alloys, 
clear short-range order is shown in agreement with available experiments for the majority 
pairs. The situation is different for the minority pairs within the Co-rich region. In this 
region, both simulated Zr-Zr and Nb-Nb RPDFs are correlated and exhibit transfers of 
atoms from the first shells of neighbours to the second and third. The Honeycutt- Andersen 
analysis exhibits mainly both distorted and pure icosahcdral orders of various degrees in 
competition with other crystalline orders in the liquid phases and supercooled alloys. 

JDG would like to thank the Government of Western Australia for a Premier's Research 
Fellowship. 



[1] A. Inouc, B. Shen, H. Koshiba, H. Kato, and A.R. Yavari. Nature Materials, 2:661, 2003. 

[2] A.R. Yavari. Nature Materials, 6:181, 2007. 

[3] D. B. Miracle. Nature Materials, 3:697, 2004. 

[4] U.K. Rofiler and H. Teichler. Phys. Rev. E, 61:394, 2000. 

[5] D. Wolf, V. Yamanakov, S.R. Phillpot, A. Mukherjee, and H. Gleitcr. Acta Materialia, 53:1, 
2004. 

[6] M.S. Daw and M.I. Baskes. Phys. Rev. B, 29:6443, 1984. 

[7] M.W. Finnis and J.E. Sinclair. Phil. Mag. A, 50:45, 1984. 

[8] Ch. Hausleitner and J. Hafner. Phys. Rev. B, 45:115, 1992. 

17 



[9] R Philips, J. Zou, A.E. Carlsson, and M. Widom. Phys. Rev. B, 49:9322, 1994. 
[10] M.W. Finnis. Interatomic Forces in Condensed Matter. Oxford University Press, Oxford, 
2003. 

[11] D.J. Oh and R.A. Johnson. J. Mater. Res., 3:471, 1988. 
[12] R.A. Johnson and D.J. Oh. J. Mater. Res., 4:1195, 1989. 
[13] R.A. Johnson. Phys. Rev. B, 39:12554, 1989. 
[14] J. Cai and Y.Y. Ye. Phys. Rev. B, 54:8398, 1996. 
[15] R.F. Zhang, Y. Kong, and B.X. Liu. Phys. Rev. B, 71:214102, 2005. 
[16] G.D. Ackland, D.J. Bacon, A.F. Calder, and T. Harry. Phil. Mag. A, 75:713, 1997. 
[17] M.I. Mendelev, S. Han, D.J. Srolovitz, G.D. Ackland, D.Y. Sun, and M. Asta. Phil. Mag., 
83:3977, 2003. 

[18] S.L. Dudarev and P.M. Derlet. J. Phys.: Condens. Matter, 17:7097, 2005. 
[19] G.D. Ackland and S.K. Reed. Phys. Rev. B, 67:174108, 2003. 
[20] G.D. Ackland. Phys. Rev. Lett., 97:015502, 2006. 

[21] G. Suran, M. Naih, M Rivoire, and J.C.S. Levy. J. Appl. Phys., 67:5649, 1990. 
[22] R.P. Gupta. Phys. Rev. B, 23:6265, 1981. 

[23] D. Tomanek, A.A. Aligia, and C.A. Balseiro. Phys. Rev. B, 32:5051, 1985. 
[24] K.W. Jacobscn, J.K. N0skkov, and M.J. Puska. Phys. Rev. B, 35:7423, 1987. 
[25] X.D. Dai, Y. Kong, J.H. Li, and B.X. Liu. J. Phys.: Condens. Matter, 18:4527, 2006. 
[26] O. Yifang, Z. Bangwei, L. Shuzhi, and J. Zhanpeng. Z. Phys. B - Condensed Matter, 101:161, 
1996. 

[27] Ch. Hausleitner and J. Hafner. J. Phys.: Condens. Matter, 2:6651, 1990. 

[28] Z. Bangwei and O. Yifang. Phys. Rev. B, 48:3022, 1993. 

[29] P. Turchi and F. Ducastelle. The Recursion Method and Rs Applications. Springer- Verlag, 
Berhn, 1985. 

[30] F. Cleri and V. Rosato. Phys. Rev. B, 48:22, 1993. 

[31] C.S Yoo, P. Soderlind, and H. Cynn. J. Phys.: Condens. Matter, 10:L311, 1998. 
[32] P. Modak, A.K. Verma, R.S. Rao, B.K. Godwal, and R. Jeanloz. Phys. Rev. B, 74:012103, 
2006. 

[33] F. Willaime and C. Massobrio. Phys. Rev. B, 43:11653, 1991. 

[34] F. Ducastelle. PhD thesis, Universite de Paris-Sud, Orsay, France, 1972. 

18 



[35] R. Pasianot, D. Farkas, and E.J. Savino. Phys. Rev. B, 43:6952, 1991. 

[36] J. Haglund, A. Fernandez- Guiller met, G Grimvall, and M. Korling. Phys. Rev. B, 48:11685, 
1993. 

[37] F. Willaime and C. Massobrio. Phys. Rev. Lett, 63:2244, 1989. 

[38] S. Baroni, A. Dal Corso, S. de Gironcoli, P. Giannozzi, C. Cavazzoni, G. Ballabio, S. Scandolo, 
G. Chiarotti, P. Focher, A. Pasquarello, K. Laasonen, A. Trave, R. Car, N. Marzari, and 
A. Kokalj. |http:// www.pwscf.orgl 

[39] J.P. Perdew, K. Burke, and M. Ernzerhof. Phys. Rev. Lett, 78:1396, 1996. 

[40] H.J. Monkhorst and J.D. Pack. Phys. Rev. B, 13:5188, 1976. 

[41] J.D. Gale and A.L. Rohl. Molecular Simulation, 29:291, 2003. 

[42] C. Kittel. Introduction to the Solid State Physics. Wiley, New York, 1966. 

[43] G. Simmons and H. Wang. Single Crystal Elastic Constants and Calculated Aggregated Prop- 
erties. MIT Press, Cambridge, 1971. 

[44] P.J. Craievich, J.M. Sanchez, R.E. Watson, and M. Weinert. Phys. Rev. B, 55(2) :787, 1997. 

[45] M.I. Baskes and R.A. Johnson. Modelling Simul. Mater. Sci. Eng., 2:147, 1994. 

[46] K. H. J. Buschow. Journal of the Less-Common Metals, 85:221, 1982. 

[47] J.K. Pargeter and W. Hume-Rothery. Journal of the Less-Common Metals, 12:366, 1967. 

[48] C.J. Smithell. Smithells Metals Reference Book. Butterworths, London, 1983. 

[49] F.B. Hildebrand. Introduction to Numerial Analysis. McGraw-Hill, 1965. 

[50] D. Holland-Moritz, T. Schenk, R. Bellissent, V. Simonet, K. Funakoshi, J.M. Merino, T. Bus- 
laps, and S. Reutzel. J. Non-Cryst Solids, 312-314:47, 2002. 

[51] G.M. Bhuiyan, M. Silbert, and M.J. Stott. Phys. Rev. B, 53:636, 1995. 

[52] X.J. Han, J.Z. Wang, M. Chen, and Z.Y. Guo. J. Phys.: Condens. Matter, 16:2565, 2004. 

[53] Y Kita, J.B. Van Zydveld, Z. Morita, and T. lida. J. Phys.: Condens. Matter, 6:811, 1994. 

[54] G. Kresse and J. Hafner. Phys. Rev. B, 48:13115, 1993. 

[55] T. Schenk, D. Holland-Moritz, V. Simonet, R. Bellissent, and D. M. Herlach. Phys. Rev. Lett., 
89:075507, 2002. 

[56] J.D. Honeycutt and H.C. Andersen. J. Phys. Chem., 91:4950, 1987. 
[57] H. Jonsson and H.C. Andersen. Phys. Rev. Lett, 60:2295, 1988. 
[58] T.H. Kim and K.F. Kelton. J. Chem. Phys., 126:054513, 2007. 
[59] N Jakse and A Pasturel. Phys. Rev. Lett, 91:195501, 2003. 

19 



[60] J.L. Rodriguez-Lopez, F. Aguilera-Granja, K. Michaelian, and A. Vega. J. of Alloys and 

Compounds, 369:93, 2004. 
[61] A.K. Shurin, P.I. Kripyakevich, and E.I. Gladyshevskii. Kristallografiya, 10(3):414, 1965. 
[62] Yu.B. Kuz'ma, V.Ya. Markiv, Yu.V. Voroshilov, and R.V. Skolozdra. Izvestiya Akademii Nauk 

SSSR, Neorganicheskie Materialy, 2:259, 1966. 
[63] A.R. Miedema, RF. de Chatel, and F.R. de Boer. Physica B, 100:1, 1980. 
[64] Yu.A. Babanov, A.F. Sidorenko, A.V. Ryazhkin, V.R. Shvetsov, J. Moessinger, and H. Kron- 

mueller. Nucl. Inst, and Meth. in Phys. Res. A, 405:400, 1998. 
[65] B.K Vainshtein, V.M. Fridkin, and V.L. Indenbom. Structure of Crystals, volume 2 of Modern 

Crystallography. Springer, 3'"'' edition, 2000. 
[66] M. Jamet, V Dupuis, P. Melinon, G. Guiraud, A. Perez, W. Wernsdorfer, A. Traverse, and 

B. Baguenard. Phys. Rev. B, 62:493, 2000. 
[67] D.R. Nelson and F. Spaepen. volume 42, page 1. Academic, Boston, 1989. 



20 




r(A) 



FIG. 1: Simulated and experimental RPDF g{r) of cobalt at a temperature of 1670 K, below the 
melting point. Experimental data from 5^] are shown by points. The insert shows the simulated 
rg{r) at 1670 K and its analysis in Gaussian peaks (dashed curves). 
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r(A) 



FIG. 2: Simulated and experimental RPDF g{r) of liquid zirconium at a temperature of 2290 K. 



Experimental data from 



55[ are shown by points. The insert shows the simulated rg{r) at 2290 K 



and its analysis in Gaussian subpeaks (dashed curves) 
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r(A) 

FIG. 3: Simulated RPDF g{r) of liquid niobium at a temperature of 2750 K. The insert shows the 
simulated rg{r) and its analysis in seven Gaussian subpeaks (dashed curves). 
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Zr(2290K) 

-- Co(1800K) 
■■■■ Nb(2750K) 




FIG. 4: Bond angle distribution at T = 2290 K for liquid Zr (solid line), at T = 1800 K for liquid 
Co (dashed line) and at T = 2750 K for liquid Nb (dotted line). The peaks in the bond angle 
distribution for perfect icosahedral order are indicated by the vertical lines. 
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FIG. 5: Radial partial distribution functions of Coo.gZro.i at 300 K compared with experiment 



(open circles) 



641 ] ■ For each atomic pair, the insert shows the simulated rg(r) at 300 K and its 



analysis in Gaussian peaks (dashed curves). 
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FIG. 6: Simulated radial partial distribution functions of Coo.gNbo.i quenched at 300K. 
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FIG. 7: Bond angular distribution functions of Coo.gZro.i and Coo.gNbo.i quenched at 300K. 
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